clear all
set more off , permanently
capture log close
graph set window fontface "Garamond"


********************************************************************************
******* This file prepares station level km to Franco-German border  ***********
********************************************************************************
	
************************************************* 
*********           GEO data           **********
************************************************* 

//	Load coordinates of Germany and France
shp2dta using "$raw\05_Coordinates\VG250_KRS.shp" ,  database("$raw\05_Coordinates\VG250_KRS_database.dta") coordinates("$raw\05_Coordinates\VG250_KRS_coord.dta") replace 
shp2dta using "$raw\05_Coordinates\VG250_STA.shp" ,  database("$raw\05_Coordinates\VG250_STA_database.dta") coordinates("$raw\05_Coordinates\VG250_STA_coord.dta") replace 

shp2dta using "$raw\05_Coordinates\FRA_adm3.shp" ,  database("$raw\05_Coordinates\FRA_adm3_database.dta") coordinates("$raw\05_Coordinates\FRA_adm3_coord.dta") replace 
shp2dta using "$raw\05_Coordinates\FRA_adm4.shp" ,  database("$raw\05_Coordinates\FRA_adm4_database.dta") coordinates("$raw\05_Coordinates\FRA_adm4_coord.dta") replace 

shp2dta using "$raw\05_Coordinates\DEU_adm3.shp" ,  database("$raw\05_Coordinates\DEU_adm3_database.dta") coordinates("$raw\05_Coordinates\DEU_adm3_coord.dta") replace 


//	Also save mercador projection
foreach i in VG250_KRS VG250_STA FRA_adm3 FRA_adm4 DEU_adm3 {
	use "$raw/05_Coordinates/`i'_coord.dta", clear
	replace _X = _X * 0.64124151 // Adjustment for map:cos(lat_mean*0.0174533); cos(50.11552*0.0174533)
	save "$dta/`i'_coord_map.dta", replace
	}

	
************************************************ 
********           Load data              ****** 
************************************************ 	

// Load France coordinates 
cd "$dta"

forval i=2013(1)2014 {
	use id_pdv lat lon using "$raw\04_France_prices\PrixCarburants_annuel_`i'.dta", clear
	bysort id_pdv: keep if _n == 1
	tempfile pc`i'
	save "`pc`i''", replace
	}
clear
forval i=2013(1)2014 {
	append using "`pc`i''"
	}
bysort id_pdv: keep if _n == 1


// Fix latitude and longitude variables
foreach i in lat lon {
	tostring `i', gen(`i'_str) force
	replace `i'_str = subinstr(`i'_str,"-.","-0.",.)
	split `i'_str, parse(".")
	
	destring `i'_str1, gen(`i'1) force
	replace `i'1 = `i'1 / 100000
	}
keep id_pdv lat1 lon1
rename lat1 latitude
rename lon1 longitude
** deal with missings later
gen longitude_map_french = longitude* 0.64124151 
save french_stations_coordinates.dta , replace
	
//	Map: Distribution french stations
use "$raw\05_Coordinates\FRA_adm3_database.dta",clear

spmap using "$dta\FRA_adm3_coord_map.dta" if NAME_1 !="Corse", ///
clmethod(unique) ///
id(_ID) legend(size(small)) ///
point(data("$dta/french_stations_coordinates.dta") ///
	select(keep if latitude > 20 & latitude <51 & (longitude < 8.5)) ///
	xcoord(longitude_map_french) ycoord(latitude) fcolor(navy)  ocolor(black) size(*0.5) legenda(off)) /* legenda(on) legtitle("") leglabel("Stations")*/
*graph export "$output/spmap_french_stations.pdf", replace



// Load Germany coordinates 
cd "$dta"
use stations.dta , clear

// Keep relevant variables
keep id_data latitude longitude
compress
save stations_coordinates.dta, replace 



//	Combined Germany & France coordinates data
use "$raw\05_Coordinates\DEU_adm3_coord.dta", clear
gen obs = _n
sort _ID obs
tostring _ID, gen(_ID_DEU)
replace _ID_DEU = "D"+ _ID_DEU
drop _ID
tempfile DEU
save "`DEU'"

use "$raw\05_Coordinates\FRA_adm3_coord.dta", clear
gen obs = _n
sort _ID obs
tostring _ID, gen(_ID_FRA)
replace _ID_FRA = "F"+ _ID_FRA
drop _ID
append using "`DEU'"

gen _ID_both =  _ID_FRA
replace _ID_both = _ID_DEU if _ID_both ==""

egen _ID = group(_ID_both)

sort _ID obs
keep _ID _X _Y 
save "$dta\FRA_DEU_adm3_coord.dta", replace
replace _X = _X * 0.64124151 // Adjustment for map:cos(lat_mean*0.0174533); cos(50.11552*0.0174533)
save "$dta\FRA_DEU_adm3_coord_map.dta", replace
	
//	Prepare coordinates data
use "$raw\05_Coordinates\VG250_STA_coord.dta", clear
gen id2 = _n
drop if _X == .
drop if _Y == .
save  "$dta\VG250_STA_coord2.dta", replace

//	French-German border coordinates data
	** north start: latitude = 49.469434, longitude = 6.367541
	** south end: latitude = 47.590036, 7.589113
	** west: longitude = 8.6
	** east:longitude 6.0
	
use "$dta\VG250_STA_coord2.dta", clear
keep if _X >= 6 & _X <= 8.6
keep if _Y >= 47.590036 & _Y <= 49.469434 
drop if _Y < 48.167643 & _X > 7.66

save "$dta\VG250_STA_coord2_frenchborder.dta", replace
replace _X = _X * 0.64124151 // Adjustment for map:cos(lat_mean*0.0174533); cos(50.11552*0.0174533)
save "$dta\VG250_STA_coord2_frenchborder_map.dta", replace

//	Find border stations
use stations_coordinates.dta , clear
append using french_stations_coordinates.dta 

gen id = _n
geonear id latitude longitude using "$dta\VG250_STA_coord2.dta", neighbors(id2 _Y  _X) ignoreself wide nearcount(1)
rename km_to_nid km_to_border
drop nid
geonear id latitude longitude using "$dta\VG250_STA_coord2_frenchborder.dta", neighbors(id2 _Y  _X) ignoreself wide nearcount(1)
rename km_to_nid km_to_french_border


//	Keep 
keep id_data id_pdv latitude longitude km_to_border km_to_french_border
order id_data id_pdv

//	Generate coordinate for map
gen longitude_map = longitude* 0.64124151 // Adjustment for map:cos(lat_mean*0.0174533); cos(50.11552*0.0174533)

//	Save
save stations_dist_to_border.dta, replace





// Gen Donut stations map: 20 to 100 km

use stations_dist_to_border.dta, clear

gen country = 0 
replace country = 1 if id_data !=""
label def country 0 "France" 1 "Germany" , replace

label val country country
geoinpoly latitude longitude using "$dta\FRA_DEU_adm3_coord.dta"

bysort _ID: gen _ID_n = _n
bysort _ID: egen km_to_border_mean = mean(km_to_border)
bysort _ID: egen km_to_french_border_mean = mean(km_to_french_border)

drop if km_to_french_border > 250
sort _ID
	
spmap country if _ID_n  == 1 & latitude>46.2 using "$dta\FRA_DEU_adm3_coord_map.dta", ///
clmethod(unique) ///
id(_ID) clnumber(8) fcolor(Blues) legend(size(small)) ///
line(data("$dta\VG250_STA_coord2_frenchborder_map.dta")  color(black) size(*3)  legenda(on) leglabel("Border") )  legend(pos(5)) ///
point(data("$dta/stations_dist_to_border.dta") xcoord(longitude_map) ycoord(latitude) fcolor(yellow)  ocolor(black) size(*1) legenda(on) legtitle("") leglabel("Stations") select(keep if km_to_french_border > 20 & km_to_french_border < 100))


graph export "$output/Figure_19_spmap_km_to_french_border_donut_20_100km.pdf", replace



// Erase intermediate data 
erase DEU_adm3_coord_map.dta 
erase FRA_adm3_coord_map.dta
erase FRA_adm4_coord_map.dta
erase FRA_DEU_adm3_coord.dta
erase FRA_DEU_adm3_coord_map.dta 
erase VG250_KRS_coord_map.dta
erase VG250_STA_coord_map.dta 
erase VG250_STA_coord2.dta 
erase VG250_STA_coord2_frenchborder.dta 
erase VG250_STA_coord2_frenchborder_map.dta 
erase stations_coordinates.dta
